Effects of phytoplankton physiology on global ocean biogeochemistry and climate

The similarity of the average ratios of nitrogen (N) and phosphorus (P) in marine dissolved inorganic and particulate organic matter, dN:P and pN:P, respectively, indicates tight links between those pools in the world ocean. Here, we analyze this linkage by varying phytoplankton N and P subsistence quotas in an optimality-based ecosystem model coupled to an Earth system model. The analysis of our ensemble of simulations discloses various feedbacks between changes in the N and P quotas, N2 fixation, and denitrification that weaken the often-hypothesized tight coupling between dN:P and pN:P. We demonstrate the importance of particulate N:C and P:C ratios for regulating dN:P on the global scale, with marine oxygen level being an important control. Our analysis provides further insight into the potential interdependence of phytoplankton physiology and global climate conditions.


INTRODUCTION
In the 1930s, A. C. Redfield found the similarity between the average ratios of nitrogen (N) and phosphorus (P) in marine dissolved inorganic and particulate organic matter (dN:P and pN:P), which suggests a tight relationship between dissolved and particulate pools (1). This similarity is difficult to explain because observations show considerable variations of pN:P and dN:P in time and space. Previous studies have focused mostly on how non-N 2 -fixing phytoplankton N:P (phyN:P) may regulate dN:P under the assumption of a globally constant phyN:P (2)(3)(4). This assumption is at odds with observations revealing that ambient conditions (nutrients, light, and temperature) strongly affect phytoplankton stoichiometry (5)(6)(7)(8). In addition, variations in phyN:P appear to contribute substantially to the pronounced spatiotemporal variability of pN:P across the world ocean (9). Thus, a constant phyN:P likely introduces considerable bias to model investigations with regard to the regulation of dN:P.
Maintaining dN:P close to the Redfield ratio could be the result of effective averaging within ecosystems on a seasonal time scale or of large-scale mixing via ocean circulation (10). Several physiological and ecological hypotheses also have been proposed to explain the similarity of dN:P and pN:P in view of the variability of phytoplankton C:N:P (11,12). Nevertheless, a convincing explanation for this phenomenon remains elusive. A useful model approach for elucidating interdependencies between dN:P and pN:P is to combine biogeochemical processes with ocean circulation, i.e., advection and mixing. From a climate perspective, changes in dN:P and pN:P affect marine primary production and the vertical export and remineralization of particulate organic matter (POM) in the ocean, which have an impact on atmospheric carbon dioxide (CO 2 ) concentration. Earth system models can be used for testing hypotheses that may explain global-scale phenomena (13), such as the relations among dN:P, pN:P, and global climate conditions. Here, we apply the optimality-based plankton ecosystem model (OPEM) that relates the optimal growth of the plankton to ambient environmental conditions while resolving variations in elemental ratios (N:C, P:C, and N:P) of the POM and dissolved inorganic nutrients (14,15). The OPEM ecosystem contains two functional types of primary producers, (non-N 2 fixing) phytoplankton and facultative diazotrophs. OPEM allocates cellular resources of primary producers for the acquisition of C, N, and P to maximize the net balanced growth rate. It is implemented in the UVic Earth system climate model (UVic-ESCM), an Earth system model of intermediate complexity that is computationally feasible for a large number of ensemble simulations (16). UVic-ESCM contains an atmospheric component and therefore allows examining mutual dependencies between marine biogeochemistry and the global climate.
The effects of variable phytoplankton C:N:P on the global biogeochemistry have received considerable attention in recent years. Variable C:N:P affects surface nutrient distributions and inventories (10,(17)(18)(19) and could also mitigate the changes in C export under different climate conditions (20)(21)(22)(23). Previous studies have considered how variations in elemental stoichiometry may be linked to the diversity of phytoplankton functional types (10,24) and how allowing for variable stoichiometry affects future climate projections (25).
Here, we revisit the more fundamental question of what drives the similarity between pN:P and dN:P on the global scale. We start out with calibrating a reference simulation for OPEM coupled to UVic-ESCM (UVic-OPEM) against data from the World Ocean Atlas 2013 (26) under preindustrial climate conditions with a fixed atmospheric partial pressure of CO 2 (pCO 2 ) (see Materials and Methods). To elucidate the linkage between pN:P and dN:P, we then perturb phyN:P by varying the subsistence N and P quotas and evaluate the resulting changes in dN:P and pN:P in a factorial sensitivity analysis. Our UVic-OPEM reference simulation describes spatial patterns of marine particulate and dissolved elemental stoichiometry, as shown earlier (14,15). The subsistence quotas largely characterize different functional types (27,28), so that different subsistence quotas can be viewed as representing the phytoplankton community by different "average" functional types, resulting in different stoichiometry patterns in the world ocean. Thus, we can investigate feedbacks between phyN:P and dN:P and the potential coupling and interdependence between pN:P and dN:P.
Our sensitivity analysis comprises an ensemble of 400 parameter sets with different combinations (20 by 20) of the N and P subsistence quotas of phytoplankton (Q N 0;phy , Q P 0;phy ). The ranges of Q N 0;phy and Q P 0;phy extend by factors between about 3.5 and 5 around the values of the reference simulation, yielding Q N 0;phy :Q P 0;phy ratios between 1 and 400 mol mol −1 . Atmospheric pCO 2 is set to evolve freely to analyze potential impacts of variations in Q N 0;phy and Q P 0;phy on the global C cycle. All other settings and parameters are left unchanged from the reference simulation, including those of the diazotrophs (table S1). Other factors that can affect phytoplankton community structure and the biogeochemical cycles of N and P, such as the source of inorganic N (29,30) and dissolved organic P (31,32), are beyond the scope of our analysis, because they are not part of OPEM. The linkage between pN:P and dN:P The subsistence quotas (Q N 0;phy and Q P 0;phy ) broadly determine the ranges of the realized phytoplankton N and P quotas (phyN:C and phyP:C) in the model. Owing to the large contribution of phytoplankton to total POM in our model (fig. S1), phyN:C, phyP:C, and phyN:P largely determine pN:C, pP:C, and pN:P (cf. Fig. 1 and figs. S2 to S4). For example, at Q P 0;phy ¼ 2 mmol mol −1 , a doubling in Q N 0;phy from 0.06 to 0.12 mol mol −1 causes median phyN:C and pN:C to rise by about 50% ( Fig. 2A and fig. S5A). Global median pN:C even doubles from 0.14 to 0.28 mol mol −1 for this change in Q N 0;phy (Fig. 2D). The sensitivity of the global median pN:C is enhanced by a strong concomitant increase in the fraction of the ocean surface area with high NO À 3 concentrations, so that the median surface NO À 3 concentration increases from ≈7 mmol m −3 for Q N 0;phy ¼ 0:06 mol mol −1 to ≈21 mmol m −3 for Q N 0;phy ¼ 0:12 mol mol −1 (Fig. 2F). Similarly, at Q N 0;phy ¼ 0:06 mol mol −1 , a doubling in Q P 0;phy from 2 to 4 mmol mol −1 causes an increase by about 50% in phyP:C ( fig. S5B) and a rise in pP:C by about 40%, depending on the ambient PO 3À 4 concentration (Fig. 2B). In contrast to the relation between Q N 0;phy and the global median pN:C, however, the effect of the same increase in Q P 0;phy on the global median pP:C appears rather subdued, amounting to only about 1% (Fig. 2C). This reduced sensitivity in global median pP:C results from a negative feedback between phyP:C and surface PO 3À 4 . That is, when Q P 0;phy increases, phyP:C and, hence, the demand for PO 3À 4 increases and phytoplankton takes up more PO 3À 4 . Thus, the sea surface area with low PO 3À 4 concentration increases, which reduces the global median pP:C (Fig. 2E).

RESULTS
The divergent responses in simulated pN:C and pP:C to changes in the subsistence quotas are associated with the different controls of dissolved inorganic PO 3À 4 and NO À 3 inventories. Because the PO 3À 4 inventory is fixed in the model, changes in model parameters only cause a spatial redistribution of PO 3À 4 within the global ocean. By contrast, the NO À 3 inventory can vary, as it is sensitive to rates of N 2 fixation and denitrification in the water column and the sediments (Fig. 2, H, M, and N), which contribute to the different sensitivities of pN:C and pP:C to changes in Q N 0;phy and Q P 0;phy , respectively.
The subsistence N and P quotas for diazotrophs (Q N 0;dia and Q P 0;dia ) are fixed in the ensemble simulations, and the ranges of Q N 0;phy and Q P 0;phy encompass Q N 0;dia and Q P 0;dia . When Q N 0;phy or Q P 0;phy is high enough, the (facultative) diazotrophs can outcompete phytoplankton even in areas with high surface NO À 3 . Diazotrophs can contribute up to 37% to total surface POM in these cases (fig. S6), although N 2 fixation always remains restricted to low NO À 3 conditions. When Q P 0;phy is high, phytoplankton N:P can be as low as 1.3 mol mol −1 , and diazotroph N:P can be as high as 53 mol mol −1 , resulting in a globally averaged pN:P of up to 28 ( fig. S7, G, H, and J).

Interdependence between dN:P and pN:P
The global dN:P is most sensitive to changes in Q N 0;phy in our ensemble simulations because of the interplay between diazotrophs and other phytoplankton. The oceanic NO À 3 inventory is driven by the balance of N 2 fixation by diazotrophs and denitrification in the water column and sediments under oxygen (O 2 )-deficient conditions. The effect of varying subsistence quotas on the NO À 3 inventory is best appreciated by combining the effect on N 2 fixation with that on denitrification, mediated by net primary production (NPP) and its effect on O 2 . Increasing Q N 0;phy raises phyN:P and the consumption of NO À 3 , which increases pN:P but lowers surface NO À 3 concentrations and, hence, promotes N 2 fixation. However, the higher N demand combined with lower surface NO À 3 eventually reduces NPP and the subsequent POM export and remineralisation at depth, which leads to less O 2 consumption and, thus, less denitrification. Both increased N 2 fixation and reduced denitrification add up, increasing the NO À 3 inventory. Decreasing Q P 0;phy increases phyN:P, NPP, and also N 2 fixation, but rising NPP and POM export also leads to higher O 2 consumption, which increases denitrification and counteracts the increase in N 2 fixation. This negative feedback, thus, limits the effect of Q P 0;phy on dN:P. While previous studies, assuming constant phyN:P, have concluded that dN:P should be linearly correlated with pN:P (3,4), this relationship turns out to be more complex in our ensemble of simulations, allowing for variable phytoplankton stoichiometry. dN:P generally does increase with higher pN:P (Fig. 3A), but we identify two notable features: (i) Simulated global average pN:P and dN:P do not form a one-to-one (injective) functional relation across our ensemble of simulated worlds, differing only in phytoplankton subsistence quotas. Instead, it allows for similar dN:P being attributable to different pN:P and vice versa (Fig. 3A), and (ii) dN:P has an upper limit of around 28 in our ensemble simulations, although pN:P can increase further. The deviant, nonlinear relationship between pN:P and dN:P can be attributed to differences in spatiotemporal variations between phyN:C and phyP:C, because they exert indirect but strong control on dN:P via O 2 (Figs. 2I and 3C). A negative feedback regulating the marine nitrogen inventory has been described previously in a qualitative manner, where O 2 acts as a biogeochemical mediator in the ocean that connects variations in pN:P to variations in dN:P (35,36). Variations in O 2 concentration reflect variations in primary production and the ensuing export and respiration, which, in turn, is influenced by changes in phyN:C and phyP:C (Fig. 2I). Once O 2 is depleted (O 2 -deficient zones), denitrification sets in, reducing the NO À 3 inventory and, hence, dN:P. In turn, changes in dN:P affect the rate of N 2 fixation of the diazotrophs, which compensates for NO À 3 losses from denitrification once a steady state is reached and a different NO À 3 inventory is established.
A key finding is that contrary to dN:P, which is driven mainly by changes in Q N 0;phy , the main driver of the variations in pN:P is Q P 0;phy . The range of realized pN:P among the simulations is narrowest when dN:P is close to the Redfield ratio (region II in Fig. 3), owing to an effective negative biogeochemical feedback between variations in PO 3À 4 and pN:P: Increasing Q P 0;phy raises the P demand for N assimilation in phytoplankton, leading to a decline in both phyN:P ( fig. S7G) and surface PO 3À 4 concentrations (Fig. 2E). This results in a competitive advantage for diazotrophs, which, owing to the high diazotroph's N:P ratio, can counteract the lower phyN:P and thereby reduce the variability of pN:P. Under these conditions, NPP is co-limited by N and P.
This negative feedback is less effective in both the low and the high dN:P ranges (regions I and III in Fig. 3), however, where Q P 0;phy exerts much stronger control on pN:P, albeit for different reasons. In the high dN:P range (region I in Fig. 3A), NPP is mostly Fe limited, and relatively high surface macronutrient concentrations (Fig. 3B) largely diminish the diazotrophs' advantage. Hence, diazotrophs have little effect on NPP and O 2 and cannot compensate for changes in phyN:P. In contrast, owing to the low NO À 3 conditions associated with the lower dN:P range (region III in Fig. 3B), variations in pN:P are mostly driven by changes in diazotroph abundance benefitting from high Q P 0;phy ( fig. S6B), also shifting the ocean from NO À 3 to PO 3À 4 limitation as Q P 0;phy increases ( Fig. 3B). As a result, pN:P, global NPP, O 2 , and NO À 3 are all highly sensitive to changes in Q P 0;phy in the lower dN:P range (region III in Fig. 3, A and C). The increase in relative diazotroph abundance with increasing Q P 0;phy , together with the low Q N 0;phy , also causes a positive relation between pN:P and Q P 0;phy in this region, which is opposite to that in region I. For dN:P close to the Redfield ratio (region II in Fig. 3), pN:P is less variable because of the negative feedback between PO 3À 4 and pP:C, which is most effective when NPP is co-limited by N and P.
Biogeochemical constraints on the changes in dN:P and pN:P Two feedbacks restrict dN:P to values below ≈28 for high Q N 0;phy , in spite of the low denitrification associated with low primary production and low particulate organic carbon (POC) export found at the upper end of the dN:P range (Fig. 4D). (i) High dN:P suppresses N 2 fixation, because the facultative diazotrophs in UVic-OPEM switch to using NO À 3 at high ambient NO À 3 concentrations. (ii) Dissolved iron (dFe) limits phytoplankton photosynthesis and N assimilation in UVic-OPEM, resulting in an inverse relation between surface dFe and NO À 3 concentrations (Fig. 2, G and H). Because increasing Q N 0;phy increases phytoplankton Fe demand, surface dFe decreases toward high Q N 0;phy , resulting in more severe Fe limitation and higher NO À 3 concentrations. Thus, Fe limitation of N 2 fixation also contributes to limiting dN:P, in spite of the low denitrification associated with the low NPP and POC export at high Q N 0;phy . Together, these constitute a negative feedback between surface NO À 3 and N 2 fixation, which limits the maximum achievable surface NO À 3 concentration.
Nitrate concentrations at the lower end of the dN:P range can be very low because of high rates of denitrification associated with the low O 2 content of the ocean under these conditions. N 2 fixation often cannot keep up with denitrification, particularly if both processes become spatially coupled. The strength of the coupling between N 2 fixation and denitrification is sensitive to the ecological niche of the diazotrophs, i.e., assumptions made with respect to the drivers of N 2 fixation (37).
Two major negative feedbacks limit the changes in pN:P at the global scale in our ensemble of simulations: (i) The negative feedback between phyP:C and surface PO 3À 4 restricts the ranges of pP:C and pN:P, and (ii) the presence of diazotrophs prevents pN:P from getting as low as that of phytoplankton at high Q P 0;phy and low Q N 0;phy . The diazotrophs also limit the increase of pN:P for high Q N 0;phy and low Q P 0;phy (fig. S7). These negative feedbacks result in pN:P staying between 6 and 38 mol mol −1 (Fig. 3A), although Q N 0;phy :Q P 0;phy ranges from 1 to 400 mol mol −1 .

Biotic controls on marine O 2 and atmospheric pCO 2
The growth of phytoplankton and diazotrophs determines the rate of CO 2 fixation, which reduces atmospheric pCO 2 due to enhanced air-sea gas flux of CO 2 (Fig. 2J). We find a strong correlation between atmospheric pCO 2 and the oceanic O 2 inventory (R 2 = 0.98; Fig. 4A). A large oceanic O 2 inventory corresponds to low remineralization and, thus, low oceanic C storage and high atmospheric pCO 2 . A small oceanic O 2 inventory corresponds to high remineralization and high oceanic C storage, leaving less CO 2 in the atmosphere. Because O 2 is the main mediator between the inventories of C and N, the strong correlation between atmospheric pCO 2 and the oceanic O 2 inventory implies a mutual dependence between atmospheric pCO 2 and the marine NO À 3 inventory. pCO 2 and NO À 3 are highly correlated (R 2 = 0.96) across our ensemble simulations (Fig. 4B). Thus, the eventual imprint of Q N 0;phy and Q P 0;phy on the NO À 3 inventory in the model translates into a concomitant, roughly proportional change in atmospheric pCO 2 , which affects the climate and ocean circulation in the model. Because surface dFe largely limits the marine C fixation, atmospheric pCO 2 and surface dFe are highly negatively correlated (R 2 = 0.98; Fig. 4C).
In the ensemble simulations, the POC export is negatively correlated with the pCO 2 (Fig. 4D), and the POC export is as high as 17 petagrams of carbon (Pg C) per year when the pCO 2 is down to 88 parts per million (ppm) and is as low as 4.6 Pg C per year when the pCO 2 is close to 400 ppm. Surface PO 3À 4 concentration has been used as an indicator for atmospheric pCO 2 because the utilization of surface PO 3À 4 is associated with C fixation (38)(39)(40). This relation becomes more complicated when considering a variable pP:C, which may amplify (38,41) or damp (23) the effect of surface PO 3À 4 changes on atmospheric pCO 2 . In the present study, the correlation between pCO 2 and surface PO 3À 4 depends on whether we are looking at variable Q N 0;phy or Q P 0;phy . For constant Q P 0;phy , the correlation is mostly positive, as both NPP and PO 3À 4 decline and, hence, both pCO 2 and surface PO 3À 4 rise with increasing Q N 0;phy . However, for constant Q N 0;phy , the correlation ranges from negative for small Q N 0;phy , associated with region III in Fig. 3A, to weakly positive for large Q N 0;phy , associated with region I in Fig. 3A (Fig. 4E). In contrast to NO À 3 (Fig. 4B), no clear relation emerges between pCO 2 and surface PO 3À 4 among the whole ensemble of simulations (Fig. 4E). Thus, varying phytoplankton subsistence quotas in UVic-OPEM break the relationship between pCO 2 and surface PO 3À 4 postulated in previous studies (38)(39)(40).
Climate constraints on the N:P of phytoplankton NPP in our ensemble is maximal for a specific range of optimal combinations of Q N 0;phy and Q P 0;phy . Subsistence quotas may be viewed as costs to phytoplankton growth, also because of the higher Fe demand when Q N 0;phy increases. Thus, it is not unexpected that the maximum globally integrated NPP occurs in the simulation with the lowest Q N 0;phy and Q P 0;phy (Fig. 2L). Nevertheless, the low pCO 2 in simulations with very low subsistence quotas lowers sea surface temperature (SST; Fig. 2K) and, hence, also NPP. For a given Q P 0;phy , the maximum NPP does not occur at the lowest Q N 0;phy , when Q P 0;phy . 0:001 mol mol À 1 (blue dashed line, Fig. 5), and for a given Q N 0;phy , the Q P 0;phy that yields the maximum NPP is not the lowest when Q N 0;phy � 0:04 mol mol À 1 (red dashed line, Fig. 5). Notably, the combination of Q N 0;phy and Q P 0;phy of our reference solution, which was calibrated against global data with a median dN:P close to the Redfield ratio, lies in the vicinity of the maximum NPP for Q N 0;phy � 0:04 mol mol À 1 (open diamond, Fig. 5).

Implication for the role of phytoplankton physiology in marine nutrient cycling and global climate
Our simulations encompass a very wide spread in pCO 2 (≈100 to >350 ppm; Fig. 2J), which, in this case, obviously originates from differences in the biological soft tissue pump (42), as phytoplankton is the only model component directly affected by the subsistence quotas. This may appear unexpected because changes in the biological pump usually are considered a minor contribution to variations in ocean C storage in the current climate change context, albeit with widely differing estimates (43,44). Analyses of anthropogenic CO 2 uptake by the ocean concern transient behavior on a time scale of decades to centuries, however, whereas our results refer to equilibrium Earth system states on a time scale of 10 4 years. In addition, wide-spread air-sea CO 2 disequilibria have recently been shown to amplify effects of a changed biological pumping to such an extent that it could be considered a major contributor to the glacial-interglacial pCO 2 difference (45).
Our model results demonstrate the importance of considering variable stoichiometry when analyzing the response of marine biogeochemistry to phytoplankton physiological variations. Therefore, the negative feedbacks stabilizing the marine nitrogen cycle in our study differ fundamentally from the "nutrient thermostat" described in earlier studies (2,17), where an assumed fixed average pN:P ratio acts as a trigger for N 2 fixation. Allowing for local deviations from the fixed pN:P ratio also can have considerable effects on dN:P (17). However, the nutrient thermostat concept would essentially imply a tight spatial association of N 2 fixation and O 2 minimum zones (46), which could result in a positive feedback between N 2 fixation and denitrification (47) and, thus, render this mechanism ineffective for stabilizing the marine nitrogen cycle (48). Our results further show that it is critical to consider also variations in C to nutrient ratios, which largely affect the growth of phytoplankton and denitrification, the major loss of nitrate in the ocean. Descriptions of the physiology of tiny plankton organisms can help elucidate profound interdependencies between particulate organic and dissolved inorganic matter pools in the global ocean. Accordingly, laboratory and field experiments that focus on physiological details, such as subsistence quotas of diverse phytoplankton, are not only valuable for explaining plankton dynamics but also of relevance for global marine biogeochemistry and climate.

Model description
For our study, we apply the UVic-ESCM (16) coupled to the OPEM (14,15), referred to as UVic-OPEM, in the following. The UVic-ESCM contains three components, a terrestrial model, a one-layer energy-moisture balance atmosphere model, and a three-dimensional ocean general circulation model. The horizontal resolution of the ocean, land, and atmosphere components is 3.6°× 1.8°, and the ocean has 19 vertical layers. UVic-ESCM is widely used as a comprehensive global C cycle model for simulations under various climate settings (49)(50)(51) and has also participated in the Climate Model Intercomparison Project 6 (52).
OPEM comprises four POM pools: non-N 2 -fixing phytoplankton (phytoplankton in the following), diazotrophs, zooplankton, and detritus. The growth of phytoplankton and diazotrophs follows an optimality assumption, where dynamic resource allocation between light-harvesting, PO 3À 4 and NO À 3 uptake, and N 2 fixation of diazotrophs maximizes growth rate under changing environmental conditions (light, temperature, and nutrient concentrations) (28).
Bioavailable Fe is explicitly resolved, and primary production can become limited by the ambient Fe concentration. This limitation is introduced as a Monod equation, with different half-saturation Fe concentrations for phytoplankton and diazotrophs.
Phytoplankton nutrient uptake and CO 2 fixation follow the temperature function in (53). Diazotrophs facultatively fix N 2 , i.e., they use NO À 3 when N 2 fixation would yield a lower net growth rate. The formulation of zooplankton foraging is the optimal current feeding model (54). The only sinking compartment is the detritus, which is fed by phytoplankton, diazotrophs, and zooplankton mortality, as well as sloppy feeding and egestion of zooplankton. The elemental composition is dynamic in phytoplankton, diazotrophs, and detritus, while the zooplankton has a fixed C:N:P ratio of 106:16:1. Details of the model description and evaluation have been reported in previous articles (14,15). For the current study, the code has been updated to UVic-ESCM version 2.10 (52), which includes the consideration of benthic denitrification (55), and we follow the original UVic temperature function for diazotrophs.

Calibration of the reference simulation
For ecosystem parameter calibration, we conducted 600 simulations under preindustrial conditions with parameter sets constructed with the Latin hypercube (LHC) method (14). The LHC sampling has been preferred over the estimation of parameter values with an iterative optimization algorithm here, mainly for two reasons: (i) All 600 simulations can be done in parallel, which reduces computational time considerably, and (ii) the applied metric can be refined and then applied to the existing ensemble without the need for extra model optimizations. Every ensemble member has been spun up for 10,000 years under preindustrial climate conditions with a prescribed atmospheric CO 2 concentration of 284.3 ppm, achieving approximate steady-state conditions. The reference parameter set (table S1) produces the lowest data model misfit, according to the minimum of a likelihood-based cost function (14). In this study, we apply a refined version of the cost function, which still considers surface chlorophyll a and PO 3À 4 , but NO À 3 and O 2 are replaced by N � ¼ NO À 3 À 16� PO 3À 4 þ 2:9 mmol m À 3 and modified apparent oxygen utilization (AOU), AOU � ¼ AOU À 2� PO 3À 4 , respectively. NO À 3 , PO 3À 4 , and O 2 data are from the World Ocean Atlas 2013 (26).

Setup of the sensitivity analysis
In OPEM, phytoplankton N and P cell quotas [phyN:C in (molN)(molC) −1 and phyP:C in (molP)(molC) −1 ] are allocated to maintenance and growth. The maintenance of a baseline physiology is associated with the subsistence quotas (Q N 0;phy and Q P 0;phy ), while the remainders (phyN:C À Q N 0;phy and phyP:C À Q P 0;phy ) are allocated to nutrient uptake, biosynthesis, and photosynthesis (28). Variations in resource allocation are the primary cause of the plasticity of phytoplankton stoichiometry in our model. Corresponding allocation factors are determined by the optimal growth strategy, and they reflect the acclimation state of phytoplankton. The subsistence N quota Q N 0;phy is proportional to the demands for structural material, and Q P 0;phy represents the P bound in the nucleus and membranes. Q N 0;phy and Q P 0;phy are treated as global constants in UVic-OPEM. Because the parameter values of the non-N 2 -fixing phytoplankton are assumed to be representative means of the entire community, any change of the subsistence quotas (Q N 0;phy and Q P 0;phy ) can be viewed as a shift in community composition.

Supplementary Materials
This PDF file includes: Figs. S1 to S7 Table S1